The field of the invention is quantitative imaging and material decomposition. More particularly, the invention relates to a method for determining the mass fractions of constituent components of an object using CT imaging and a post-reconstruction material-basis model.
In a computed tomography system, an x-ray source projects a fan-shaped beam which is collimated to lie within an X-Y plane of a Cartesian coordinate system, termed the “imaging plane.” The x-ray beam passes through the object being imaged, such as a medical patient, and impinges upon an array of radiation detectors. The intensity of the transmitted radiation is dependent upon the attenuation of the x-ray beam by the object and each detector produces a separate electrical signal that is a measurement of the beam attenuation. The attenuation measurements from all the detectors are acquired separately to produce the transmission profile at a particular view angle.
The source and detector array in a conventional CT system are rotated on a gantry within the imaging plane and around the object so that the angle at which the x-ray beam intersects the object constantly changes. A group of x-ray attenuation measurements from the detector array at a given angle is referred to as a “view”, and a “scan” of the object comprises a set of views acquired at different angular orientations during one revolution of the x-ray source and detector. In a 2D scan, data is processed to construct an image that corresponds to a two dimensional slice taken through the object. The prevailing method for reconstructing an image from 2D data is referred to in the art as the filtered backprojection technique. This process converts the attenuation measurements from a scan into integers called “CT numbers” or “Hounsfield units”, which are used to control the brightness of a corresponding pixel on a display.
Dual source CT systems have two separate x-ray sources and associated detector arrays, which rotate together in the gantry during a scan. The x-ray sources may be operated at different energy levels to acquire two image data sets from which a low energy and a high energy image may be reconstructed.
Quantitative imaging using CT systems has experienced tremendous growth in recent years, in terms of both the basic technology and new clinical applications. CT-based quantitative imaging exploits differences in x-ray attenuation between different materials. In CT images, because different materials cause different levels of x-ray scattering and absorption, proper calibration of image pixel values versus x-ray beam energy can be used to qualitatively and quantitatively evaluate an imaged object's material composition.
The degree to which a given material blocks x-ray transmission, is measured by an attenuation coefficient, which accounts for both energy absorption and the scattering of photons. Often, mass attenuation coefficients, which measure attenuation per unit mass, are utilized, because they do not change with the density of the material. Tabulated mass attenuation coefficients for different elements (Z=1˜92) are readily available in the database of the National Institute of Standards and Technology (NIST). The NIST database also includes 48 compounds and mixtures, covering nearly all tissues found in the human body. The mass attenuation coefficient for compounds and mixtures having than more than one material is simply the summation of weighted mass attenuation coefficients of each constituent material, wherein the weighting factor is the mass fraction of each constituent material. Material decomposition techniques can be employed to calculate these mass fractions using known mass attenuation coefficients and dual energy CT measurements. In principle, this can only be done for objects having two constituent materials, as dual energy CT provides only two independent measurements.
Alternatively, material decomposition techniques can quantify and object under investigation by analyzing the physical mechanisms that cause attenuation. For the x-ray energies in the medical diagnostic range, the mechanisms responsible for material attenuation are the photoelectric effect and Compton scattering, which can be approximately modeled using effective atomic number, density, and x-ray energy information. Therefore, instead of obtaining the mass fraction of each material, the effective atomic number and density of the imaged object can be determined using a model of these two mechanisms and dual energy CT measurements.
In 1976, Alvarez and Macovski proposed a method to couple the attenuation coefficient model with CT measurements in order to determine the atomic number and density of a material. First, the attenuation coefficient (μ( . . . )) is modeled as a linear combination of the photoelectric effect and Compton scattering, as follows:
                              μ          ⁡                      (                          x              ,              y              ,              E                        )                          =                                                            a                1                            ⁡                              (                                  x                  ,                  y                                )                                      ⁢                          1                              E                3                                              +                                                    a                2                            ⁡                              (                                  x                  ,                  y                                )                                      ⁢                                                            f                  KN                                ⁡                                  (                  E                  )                                            .                                                          Eqn        .                                  ⁢        1            
In Eqn. 1, the photoelectric effect is inversely proportional to the energy level (E) cubed, and the Compton scattering is modeled by Klein-Nishina formula. The terms a1(x, y) and a2(x, y) are related to the atomic number and physical density of the materials under investigation. For a CT scan, this model is expressed as the following line integral:
                                                        ∫                                                μ                  ⁡                                      (                                          x                      ,                      y                      ,                      E                                        )                                                  ⁢                                  ⅆ                  s                                                      =                                                            A                  1                                ⁢                                  1                                      E                    3                                                              +                                                A                  2                                ⁢                                                      f                    KN                                    ⁡                                      (                    E                    )                                                                                ;                ⁢                                  ⁢                  where          ⁢                      :                                              Eqn        .                                  ⁢        2                                                                    A              1                        =                          ∫                                                                    a                    1                                    ⁡                                      (                                          x                      ,                      y                                        )                                                  ⁢                                  ⅆ                  s                                                              ;                ⁢                                  ⁢        and                            Eqn        .                                  ⁢        3                                          A          2                =                  ∫                                                    a                2                            ⁡                              (                                  x                  ,                  y                                )                                      ⁢                                          ⅆ                s                            .                                                          Eqn        .                                  ⁢        4            
From Eqn. 2, at least two equations are needed to obtain the single solution of the two unknowns. Because dual-energy CT images the object at two different energy levels, it satisfies this requirement and the equation can be written as follows:
                              I          i                =                  ∫                                                    S                i                            ⁡                              (                E                )                                      ⁢                          ⅇ                                                -                                                            A                      1                                                              E                      3                                                                      -                                                      A                    2                                    ⁢                                                            f                      KN                                        ⁡                                          (                      E                      )                                                                                            ⁢                                          ⅆ                                  E                  ⁢                                                                          (                                                            i                      =                      1                                        ,                    2                                    )                                            .                                                          Eqn        .                                  ⁢        5            
However, it is very difficult to solve Eqn. 5 for A1 and A2. Alvarez and Macovski therefore used the following power series to approximate the integral equation:ln I1=b0+b1A1+b2A2+b3A12+b4A22+b5A1A2+b6A13+b7A23  Eqn. 6;andln I2=c0+c1A1+c2A2+c3A12+c4A22+c5A1A2+c6A13+c7A23  Eqn. 7.
The sets of coefficients {bi} and {ci} are determined by calibrations. Once the A1 and A2 are solved using Eqn. 6 and 7, a1(x,y) and a2(x,y) can be reconstructed by one of the reconstruction methods, such as filtered back projection.
This method, typically referred to as the basis-spectral method, was the first theoretical analysis on material-selective imaging using dual-energy CT. The drawback of the basis-spectral method is that it is not very accurate due to the intrinsic difficulty in modeling the photoelectric effect and Compton scattering, especially for discontinuous absorption edges. Although the basis-spectral method accounts for the photoelectric effect and Compton scattering by creating a photoelectric effect and Compton scattering map, this type of technique is more often used to give the object's effective atomic number (Z) and density (ρ).
In 1986, Kalender et al. proposed that any material's mass attenuation coefficient can be expressed as a linear combination of the coefficients of two so-called basis materials, as follows:
                                          (                          μ              ρ                        )                    ⁢                      (            E            )                          =                                                                              a                  1                                ⁡                                  (                                      μ                    ρ                                    )                                            1                        ⁢                          (              E              )                                +                                                                      a                  2                                ⁡                                  (                                      μ                    ρ                                    )                                            2                        ⁢                                          (                E                )                            .                                                          Eqn        .                                  ⁢        8            
For a CT measurement, this is expressed using the following line integral:
                                                        ∫                                                μ                  ⁡                                      (                                          x                      ,                      y                      ,                      E                                        )                                                  ⁢                                  ⅆ                  s                                                      =                                                                                                      A                      1                                        ⁡                                          (                                              μ                        ρ                                            )                                                        1                                ⁢                                  (                  E                  )                                            +                                                                                          A                      2                                        ⁡                                          (                                              μ                        ρ                                            )                                                        2                                ⁢                                  (                  E                  )                                                              ;                ⁢                                  ⁢        where                            Eqn        .                                  ⁢        9                                                                    A              1                        =                          ∫                                                                    ρ                    1                                    ⁡                                      (                                          x                      ,                      y                                        )                                                  ⁢                                  ⅆ                  s                                                              ;                ⁢                                  ⁢        and                            Eqn        .                                  ⁢        10                                          A          2                =                  ∫                                                    ρ                2                            ⁡                              (                                  x                  ,                  y                                )                                      ⁢                                          ⅆ                s                            .                                                          Eqn        .                                  ⁢        11            
This method is called the basis-material method. Similar to the basis-spectral method, dual-energy CT measurements are needed to solve the two unknowns A1 and A2. The assumption used with the basis-material method is that the attenuation coefficients of the two basis materials are known. From this assumption and the dual-energy CT measurements, the line integral equation can be written as follows:
                              I          i                =                  ∫                                                    S                i                            ⁡                              (                E                )                                      ⁢                          ⅇ                              [                                                                            -                                                                                                    A                            1                                                    ⁡                                                      (                                                          μ                              ρ                                                        )                                                                          1                                                              ⁢                                          (                      E                      )                                                        -                                                                                                              A                          2                                                ⁡                                                  (                                                      μ                            ρ                                                    )                                                                    2                                        ⁢                                          (                      E                      )                                                                      ]                                      ⁢                                          ⅆ                E                            .                                                          Eqn        .                                  ⁢        12            
Instead of solving Eqn. 12 directly, the basis-material method uses a table lookup procedure to solve the equation and produce an output that can be interpreted as components in a two-dimensional vector space, wherein the basis materials define the basis vectors. As a result, the basis-material method using dual energy CT is more efficient and clinically practical than the above-described basis-spectral approach.
Although the basis-spectral and basis-material methods differ in their modeling of the attenuation coefficients, they both belong to the “pre-reconstruction” class of methods. That is, both methods are performed with “raw data,” prior to image reconstruction. In contrast, a post-reconstruction method would be capable of analyzing reconstructed images directly.
In 2003, Heismann et al. proposed a general post-reconstruction method for performing material decomposition using CT. Under the Heismann method, the effective attenuation coefficient is first defined as follows:
                                             μ            eff                    =                    ⁢                                                    lim                                  d                  ->                  0                                            ⁢                              [                                                      -                                          1                      d                                                        ⁢                                      ln                    ⁡                                          (                                              I                                                  I                          0                                                                    )                                                                      ]                                      =                                          lim                                  d                  ->                  0                                            ⁢                                                [                                                            -                                              1                        d                                                              ⁢                                          ln                      (                                                                        ∫                                                                                    S                              ⁡                                                              (                                E                                )                                                                                      ⁢                                                          D                              ⁡                                                              (                                E                                )                                                                                      ⁢                                                          ⅇ                                                                                                -                                                                      μ                                    ⁡                                                                          (                                      E                                      )                                                                                                                                      ⁢                                d                                                                                      ⁢                                                          ⅆ                              E                                                                                                                                ∫                                                                                    S                              ⁡                                                              (                                E                                )                                                                                      ⁢                                                          D                              ⁡                                                              (                                E                                )                                                                                      ⁢                                                          ⅆ                              E                                                                                                                          )                                                        ]                                .                                                                          Eqn          .                                          ⁢          13                    
Then, as follows:
                                                        μ              eff                        =                          ∫                                                w                  ⁡                                      (                    E                    )                                                  ⁢                                  μ                  ⁡                                      (                    E                    )                                                  ⁢                                  ⅆ                  E                                                              ;                ⁢                                  ⁢        where                            Eqn        .                                  ⁢        14                                          w          ⁡                      (            E            )                          =                                                            S                ⁡                                  (                  E                  )                                            ⁢                              D                ⁡                                  (                  E                  )                                                                    ∫                                                S                  ⁡                                      (                    E                    )                                                  ⁢                                  D                  ⁡                                      (                    E                    )                                                  ⁢                                  ⅆ                  E                                                              .                                    Eqn        .                                  ⁢        15            
The effective attenuation coefficient (μeff) is determined from the CT image data and S(E) and D(E) are the tube spectrum and detector sensitivity, respectively. Like the basis-spectral method, the Heismann method treats μ(E) as the following linear combination of the photoelectric effect and Compton scattering:
                              (                                                                      μ                                      eff                    ⁢                                                                                  ⁢                    1                                                                                                                        μ                                      eff                    ⁢                                                                                  ⁢                    2                                                                                )                =                              ρ            ⁡                          (                                                                                          ∫                                                                                                    w                            1                                                    ⁡                                                      (                            E                            )                                                                          ⁢                                                  (                                                                                                                    μ                                photo                                                            ρ                                                        +                                                                                          μ                                Compton                                                            ρ                                                                                )                                                ⁢                                                  ⅆ                          E                                                                                                                                                                                ∫                                                                                                    w                            2                                                    ⁡                                                      (                            E                            )                                                                          ⁢                                                  (                                                                                                                    μ                                photo                                                            ρ                                                        +                                                                                          μ                                Compton                                                            ρ                                                                                )                                                ⁢                                                  ⅆ                          E                                                                                                                                )                                .                                    Eqn        .                                  ⁢        16            
Again like the basis-spectral method, the Heismann method models the photoelectric effect and Compton scattering as functions of atomic number and x-ray energy, respectively, as follows:
                                                                        μ                photo                            ρ                        =                          α              ⁢                                                Z                  3                                                  E                  3                                                              ;                ⁢                                  ⁢        and                            Eqn        .                                  ⁢        17                                                                    μ              Compton                        ρ                    =          β                ;                            Eqn        .                                  ⁢        18            
where α and β are constants. Therefore, Heismann's post-reconstruction method determines an object's effective atomic number and density directly from CT images. The Heismann method, though able to utilize CT data post-reconstruction, shares many of the assumptions of the above-described pre-reconstruction spectral-basis method. As a result, the Heismann method suffers from similar drawbacks.
Heismann did not determine how to use the data in both the pre-reconstruction space and post-reconstruction space, how to solve the linear equations more accurately, and other measures, such as the beam hardening effect, necessary to successfully implement such a method.
In addition, the above-described methods are incapable of performing three material decomposition using dual energy CT data. This can be problematic, as clinicians often encounter situations in which more than two materials coexist in an object, for example, when imaging bone, tissue, and an iodinated contrast material. In such cases, additional information is required because n>m. Therefore, three equations are theoretically needed to solve for three unknowns and perform three-material decomposition. Accordingly, one method for performing three-material decomposition includes employing triple-energy CT systems so that m=n. However, the mechanisms of photon attenuation are mainly the photon-electric effect and Compton scattering, triple-energy CT systems are not very useful.
A method for performing three-material decomposition using dual-energy CT data was proposed previously. The previous method assumes that the volume fraction of three materials in an imaged object can be expressed as follows:f1+f2+f3=1  Eqn. 19;
where f1 to f3 are the volume fractions as illustrated in FIG. 1, which depicts a three material object including a first material volume fraction f1 having a known density ρ1, a second material volume fraction f2 having a known density ρ2, and a third material volume fraction f3 having a known density ρ3. Therefore the effective density ρ of the three material mixture can be expressed as follows:ρ=f1ρ1+f2ρ2+(1−f1−f2)ρ3  Eqn. 20.
Thus, the effective attenuation coefficient is as follows:μ=f1ρ1+f2μ2+(1−f1−f2)μ3  Eqn. 21;
where μ1, μ2 and μ3 are the attenuation coefficients of the three materials, respectively. Accordingly, only two unknowns, f1 and f2, need to be solved to perform three material decomposition using dual energy CT data.
The volume-conservation three-material decomposition method works well in some cases, especially when the three materials are solid with clear boundaries between each other. However, there are many clinical situations where the three materials are not completely solid and lack clearly defined boundaries. The volume-conservation method is not applicable under these situations.
Accordingly, it would be desirable to have a system and method for performing three-material decomposition without the need for three-energy measurements. In addition, this system and method should provide more accurate and widely applicable clinical data than the volume-conservation method.